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Abstract 

Specific forms for both the Gibb’s and complementary dissipation potentials are chosen such 
that a complete (i.e., fully associative) potential based multiaxial unified viscoplastic model is ob- 
tained. This model possess one tensorial internal state variable, that is associated with dislocation 
substructure, with an evolutionary law that has nonlinear kinematic hardening and both thermal 
and strain induced recovery mechanisms. A unique aspect of the present model is the inclusion of 
non-linear hardening through the use of a compliance operator, derived from the Gibb’s potential, 
in the evolution law for the back stress. This non-linear tensorial operator is significant in that it 
allows both the flow and evolutionary laws to be fully associative (and therefore easily integrated) 
and greatly influences the multiaxial response under non-proportional loading paths. In addition 
to this nonlinear compliance operator, a new consistent, potential preserving, internal strain un- 
loading criterion has been introduced to prevent abnormalities in the predicted stress-strain curves, 
which are present with nonlinear hardening formulations, during unloading and reversed loading of 
the external variables. Specification of an experimental program for the complete determination of 
the material functions and parameters for characterizing a metallic matrix, e.g., TIMETAL 2 IS, is 
given. The experiments utilized are tensile, creep, and step creep tests. Finally, a comparison of this 
model and a commonly used Bodner-Partom model is made on the basis of predictive accuracy and 
numerical efficiency. 

Keywords: viscoplasticity, nonlinear hardening, isothermal, deformation, multiaxial, correlations, 
predictions 
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1 Nomenclature 


Invariants 

fi 

$ 

F 

G,y 

h 

h 

Hi 


complementary dissipation potential 

Gibb’s complementary potential 

Bingham-Prager threshold function 

normalized second invariant functions 

second invariant of effective deviatoric stress tensor 

second invariant of internal deviatoric stress tensor 

invariant material functions 


Stresses 

<rij 

Si: 

Si; 

a 7 

a «'j 

a V 

K } K 0 

Y 


Cauchy stress tensor 
deviatoric stress tensor 
effective deviatoric stress tensor 
internal state variables (stress-like) 
internal (or back) stress tensor 
deviatoric internal stress tensor 
drag stress and initial drag stress 
yield stress 


Strains , 

P .. € t total, inelastic, reversable. and thermal strain tensors, respectively 

conjugate internal state variables (displacement-like) 

^4^. internal strain tensor 


Material Parameters 

CijhjEijki 

Bq^Bi^Rcc 

P 

n,p,q 
9(G) J(F) 

H(a y ),Z(T) 

t),E,C 


elastic compliance and stiffness tensor, respectively 
hardening and thermal recovery material parameters 
material parameters 

denotes the extent of strain induced recovery 
material exponents 
material functions 
material functions 

coefficient of thermal expansion, Young’s modulus, shear modulus 


Miscellaneous 


Tref 
sgn( ) 

Qijkl, Lijkl 
hi 
0 

Hv[ ] 

(• ) 

S 

G,i> 


reference temperature 

sign function: takes on 1 or -1 depending upon the sign of argument 

internal compliance and internal stiffness operators, respectively 

Kronecker delta function 

Macauley bracket 

Heaviside unit function 

time derivative (or rate) notation 

instantaneous state trajectory 

components of the state trajectory, S • 
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2 Introduction 


A number of Titanium Matrix Composite (TMC) systems are currently being researched and evaluated 
for high temperature air frame and propulsion system applications. As a result, numerous computational 
methodologies for predicting both deformation and life are under development. An integral part of these 
methodologies is an accurate and computationally efficient constitutive model for the metallic matrix 
constituent. Furthermore, because of the proposed elevated operation temperatures for which these 
systems are designed, the required constitutive models must account for both time-dependent and time- 
independent deformations. To accomplish this we will employ a recently developed complete potential 
based framework[l] utilizing internal state variables which was put forth for the derivation of reversible 
and irreversible constitutive equations. This framework, and consequently the resulting constitutive 
model, is termed complete because the existence of the total (integrated) form of the Gibbs comple- 
mentary free energy and complementary dissipation potentials are assumed a priori. The specific forms 
selected here for both the Gibbs and complementary dissipation potentials result in a fully associative, 
multiaxial, isothermal, unified viscoplastic model with nonlinear kinematic hardening. Thus this model 
constitutes one of many models in the GVIPS (Generalized Viscoplasticity with Potential Structure) 
class of inelastic constitutive equations which can be constructed using the generalized framework of 
Arnold and Saleeb[l]. 

The particular unified GVIPS model of interest in this study possesses one tensorial internal state 
variable (i.e., the back or internal stress) that is associated with dislocation substructure and an evo- 
lutionary law that has nonlinear kinematic hardening and both thermal and strain induced recovery 
mechanisms. A unique aspect of the present model is the inclusion of non-linear hardening through the 
use of a compliance operator (derived from the Gibb’s potential) in the evolution law for the back stress. 
This non-linear tensorial operator is significant in that it allows both the flow and evolutionary laws to 
be fully associative (and therefore easily integrated) [2] and greatly influences the multiaxial response 
under non-proportional loading paths [1],[3],[4]. In addition to this nonlinear compliance operator, a 
new consistent, potential preserving, internal unloading criterion has been introduced to prevent abnor- 
malities in the predicted stress-strain curves, which are present with nonlinear hardening formulations, 
during unloading and reversed loading. This new criterion is motivated by considering the trajectory of 
the state point through the state space during unloading and hinges upon the introduction of internal 
state discontinuity inequalities for the conjugate back (or internal ) strain rate. 

A primary objective of the present study is to specify material functions and characterize the associ- 
ated material parameters for the current kinematic, isothermal GVIPS model given TIMETAL 2 IS 1 , an 
advanced titanium-based matrix commonly used in TMCs. Although both long and short term behavior 
is important, capturing the short term (or transient) behavior and rate sensitivity of the material is of 
primary importance given that the applications of interest are primarily those involving processing and 
propulsion systems. Furthermore, to illustrate the effectiveness of the present model a comparison of 
the GVIPS model and a particular Bodner-Partom (BP) model discussed in [5] (of the commonly used 
BP models [6]-[8] in the literature) is made on the basis of predictive accuracy and numerical efficiency. 

The paper begins by briefly summarizing the complete potential structure, followed by a multiaxial 
and uniaxial statement of the newly proposed GVIPS model, including a discussion on the recent addition 
of a potential preserving internal unloading criterion. A discussion regarding the characterization of the 
proposed model is then followed by numerous results illustrating the predictive accuracy and numerical 
efficiency of the model. 

3 Complete Potential Structure 

Here the basic thermodynamic framework put forth by Arnold and Saleeb [1] is summarized. Expressions 
for the Gibb’s thermodynamic and the complementary dissipation potential functions are assumed in 

1 TIMETAL 2lS is a registered trademark of TIMET, Titanium Metals Corporation. 
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terms of a number of state and internal variables characterizing the changing internal structure of the 
material. For conciseness, the discussion is limited to a case involving i) small deformations (in which 
the initial state is assumed to be stress free throughout), ii) an initially isotropic material, iii) isothermal 
conditions, and iv) the specialized potential framework discussed in [1]. A Cartesian coordinate reference 
frame and index notation are utilized (repeated Roman subscripts imply summation). 

Given the Gibb’s potential in the following form 


$ = $(a I j,Q 7 ,T,ef i ) (1) 

and assuming a priori that the inelastic strain is an independent parameter (and not an internal state 
variable), for example 


$ = E(crij) - H (o-r) - Z ( T ) - ~Y T] {T - To), 

(2) 

an expression for the total strain rate can be obtained by differentiating, that is, 


d ,-d$. _ . .i.glrr 

e,; ” dt ( daj ) ~ c,jrj<Trj + 11 +6 ' } z T 

(3) 

as well as the rate of change of the conjugate internal variables (j4(), 


• d.-d _ . 

A( - 7i ( dc ( > 

(4) 

”' h “ e „ _ -9.* . 

d<Tijd<T r , d<Tijd<T rs 

(5) 

-d'Q - d 2 H (a 7 ) 

da {da-f da^da^ 

(6) 

are the external and internal compliance operators, respectively. Note the three terms in equation (3) 
may then be identified as the reversible, irreversible (inelastic), and thermal expansion components of 
the total strain rate, respectively. Thus, 

p . . — -L. } f , -L 

(7) 

where 


f/y = Cijrs&rs 

(8) 

and 


r>\. 

II 

c* 

COl-3 

(9) 

but ifj (the inelastic strain rate) is defined separately using the concept of a 
Ci(<rij,<Xy,T). 

Given 

Q = Q(a { j ,a y ,T) 

dissipation potential 
(10) 

and using the Clausius-Duhem inequality; the flow law becomes 


,r dQ 

£,} ~ dan 

(11) 


and the evolutionary laws for the thermodynamic conjugate internal state variables: 
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( 12 ) 



Utilizing equation (4) the internal constitutive rate equations for the internal state variables is obtained, 

a 7 = As 

where 

= Ml *!" 1 = 

Thus, equations (11) and (12) represent the flow and evolutionary laws, for an assumed Q = fi(<r,j,a 7 ,T), 
and equation (13) the internal constitutive rate equations, given a Gibb’s potential $, wherein both 
potentials are directly linked through the internal state variables a 7 . 

4 Viscoplastic Constitutive Model 

A complete multiaxial statement of a GVIPS model can be derived by using the above framework, given 
a specific form for both the Gibb’s potential, and the complementary dissipation potential, 0. Form 
invariance (objectivity) of these potentials requires that they depend only on certain invariants of their 
respective tensorial arguments (i.e., an integrity basis [9] ). In the spirit of von Mises and because of the 
deviatoric nature of inelastic deformation, only the quadratic invariant will be considered at this time 
in specifying the dissipation potential. Similarly, only the linear elastic strain energy contribution will 
be considered in specifying the Gibb’s potential, with the internal state groupings being functions of the 
respective quadratic invariants. Finally, although equation (10) indicates that an unlimited number of 
internal state variables can be specified, here our attention will be restricted to a GVIPS model with 
a single independent , evolving, internal state variable. This internal state variable, a, j, is taken to be 
a second order symmetric traceless tensor that represents the internal (or back) stress associated with 
dislocation substructure. Two additional scalar state variables, one representing the drag strength (k) 
and taken to be non-evolving and the other a yield stress (Y) that is uniquely assumed to implicitly 
evolve with internal stress are considered. Both variables are associated with dislocation density and are 
included in anticipation of the subsequent nonisothermal GVIPS model. However, in the forthcoming 
nonisothermal formulation k will explicitly be taken as an evolving internal state variable with its own 
associative evolution law (see [3]). 

Consequently, the Gibb’s potential may be written as 


$ = " Crskl<7rs<Tkl 

- v i3 e! 3 + Hi(G) - Z(T) - - T 0 ) 

(15) 

and the dissipation potential as 

' to 

II 

f f(F)dF+R a B 0 J g(G)dG 

(16) 

where 


II 

1 

(17) 


Y = (l -gs/G) 

(18) 



(19) 


(13) 

(14) 
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Hi(G) = -Bq{G + B 1 CP) 


( 20 ) 


r 3 
h — 2 a ij a ij 

(21) 

J 2 — gEijSji 

(22) 

Sij = Sij — Ujj 

(23) 

Sij — <T%j — — &kk&ij 

(24) 

a,*j = a*j — 

(25) 


Note that in the preceding expression for the dissipation potential, the stress dependence, both 
external and internal, enters through the scalar functions F and G in the form of effective (£*j) and 
internal (a^) deviatoric stresses, respectively. Furthermore, the function F acts like a threshold surface, 
because when F < 0, no inelastic strain can occur. Clearly, this threshold value is dictated by the 
magnitude of both the drag strength (/c) and yield stress (Y). A unique aspect of this model is that the 
internal variable representing the yield stress is specifically taken as a special scaled function of the back 
stress and drag strength. Consequently, this allows for 1) the model to possess features of a model with 
three internal variables yet without any additional computational cost, and 2) the presence of an induced 
strain recovery term (as opposed to the common 5 ad-hoc’ introduction of such terms) in the evolution of 
the back strain (i.e., the associated conjugate variable, *4 7 ). It is important to realize that the product 
kY constitutes the radius of the initial threshold surface (see Fig la), thereby dictating from physical 
arguments, that both Y and k be always positive valued. Furthermore, given the specified form of Y in 
equation (18), it is clear that the material parameter ratio (/?/*) w iU dictate the limit value (i.e., when 
Y— 0) of the internal stress (a i; ), or the cut-off limit for dynamic recovery. 

By selecting the preceding scalar functions, a general yet complete potential-based model, with asso- 
ciated flow and evolutionary laws can be constructed. The second invariants , J 2 and h , are also scaled 
for tension. These invariants could just as easily have been scaled for shear by replacing the coefficient 
3/2 with 1/2, and modifying the definition of the magnitude of the inelastic strain rate that follows. 
Also, it should be stated that the linear form of F and Y in equations (17) and (18) was chosen in order 
to allow algebraic manipulation and analytical solution of the resulting expressions (e.g., inversion of the 
flow law), so as to ease the characterization stage of the model, as discussed subsequently in sections 4.1 
and 4.2. Note, however, that the quadratic form of F and Y does provide some additional interesting 
features in the model as discussed previously by Arnold et al. [3]. 

Finally, it is interesting to note the close similarity between the present assumed form of the threshold 
surface F of equation (14) to that used by Lemaitre and Chaboche [13] in the context of inviscid plasticity. 
The significant difference between the two forms is that Lemaitre and Chaboche [13] utilized a non- 
associated format; the form similar to equation (14) was taken to affect only the directions of inelastic 
and internal state, whereas a classical J 2 form was employed to control yielding (i.e., the consistency 
condition). In the present fully associative case, the same threshold surface F is used to affect both 
yielding, and the direction of inelastic and internal state. 

Taking the appropriate derivatives of both and fi as indicated in equations (2) through (14), one 
obtains the multiaxial isothermal specification particular to the present constitutive model. Here the 
decomposition of the total strain rate is that of equation (7), where the reversible strain rate is given by 
equations (8) and (9) and the irreversible (or inelastic) component is defined by the following flow law: 

if F -° (26) 
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Figure 1: Geometric interpretation of the internal unloading criterion relative to positions of the thresh- 
old surfaces. 


or 

=0 if F < 0 

where 


IK6II = vfii = ^ 

The internal constitutive rate equation is always given by 


(27) 


(28) 


aij - Lijr, At, (29) 

while the evolutionary law for the back strain rate is given by : 

Mi = & ~ ~ if ajjEjj > 0 (30) 

during internal loading and 

At, = Qr, lm E tmnp (ej, - lt\\i\.\\^LHv\Y) - if < 0 (31) 

during internal unloading, as indicated by the inequality a «i Ey < 0 * The internal stiffness operator is 
defined as 


Lijr s — [Qijrs] — 


3Bo(l + B\pG p 


-1) ( 7, ' jrJ K 


3Si(p-l)GP 


-2 


2(1 + B iP Gp-'(6 P - 5)) 




(32) 
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with Iijrs — fiir&js ■ 

Equation 31 constitutes a new consistent, potential preserving, internal unloading criterion that 
will prevent the classical abnormalities [13] in the cyclic response associated with nonlinear hardening 
formulations, upon unloading and reversed loading of the external variables. This new criterion is 
motivated by considering the trajectory of the state point within the state space during unloading and 
hinges upon the introduction of an internal state discontinuity inequality for the conjugate internal (or 
back) strain rate. Mathematically, this is accomplished by the introduction of a dissipation potential that 
is distinctly different when the material is undergoing internal loading versus unloading. The geometric 
interpretation of this criterion (discontinuity inequality) in the tt— plane is illustrated in Fig. 1, where 
the situation of the deformed material a tJ ^ 0, is depicted. Note that the projection of the applied 
stress vector <j ij on the tt— plane is, of course, Sij , the deviatoric stress and that the projection of 
the cylindrical surfaces of constant inelastic strain rate (F=constant) appear as concentric circles in this 
plane, although only one is illustrated in Fig. 1. The current effective stress vector = Sij — 
emanates from the translated center of the nested circles and the corresponding inelastic strain rate 
vector ifj , is collinear with it. Now, if the applied stress is reduced to a value less than the internal 
stress, the effective stress will change sign and unloading internal to the material will take place, thus 
resulting in a readjustment (state recovery, i.e., shift toward the origin of the w- plane) in internal stress 
dij as shown in Fig. lb. If the stress reversal were continued, the material would also continue to rapidly 
readjust until once again internal reloading commenced, as depicted in Fig. lc. 

Justification for this criterion stems from the work of Orowan [14] and others, and has indicated 
that upon stress reversals dislocations are remobilized and consequently rapid rearrangement of these 
dislocations is possible within the wake of the previous load path. To describe this rapid motion of 
dislocation remobilization (readjustment of internal stTess) during material unloading we introduce, in 
the spirit of Robinson [15] and Onat [11] , distinct regions within the state space in which the rate of 
the conjugate internal state variable (internal strain) changes discon tinuously. Currently, until more 
exploratory tests can be performed to fully describe the affected regions of the state space, the internal 
(back) strain rate during internal unloading is taken to be proportional to the back strain rate during 
loading, through the product of the external stiffness and internal compliance operators, see eq. (30). 
Thus, implying that the internal stress and strain during internal unloading are related by the external 
stiffness tensor; the accuracy' of such an assumption can be verified by conducting a number of transient 
stress dip tests [10]. Furthermore, an examination of the evolution of the conjugate of the internal 
stress (i.e., the internal strain, Aij), clearly shows that equation (30) possesses, as typically assumed 
in the literature [17], [18] competitive mechanisms consisting of a hardening term (which accounts for 
strengthening mechanisms) and two state recovery terms (which account for softening mechanisms). 
The first state recovery term evolves with inelasticity; it is strain induced, and is commonly called 
dynamic recovery; whereas, the second term, interchangeably called static or thermal recovery, evolves 
with time and is thermally induced. Again, it is important to note that as the internal stress evolves 
and finally reaches a limit value (dictated by the material parameters /? and /c) the yield stress will 
vanish, as will the strain induced recovery term. The added fl exibilit y provided by the dynamic recovery 
term is particularly advantageous under conditions where thermal recovery is relatively small (i.e., low 
homologous temperatures). 

Finally, the above expressions are further specialized by assuming the nonlinearity of inelasticity and 
the thermal recovery process to be represented by power law functions, that is: 

f(F) = F n (33) 

g(G) = G q (34) 

4.1 Uniaxial Simplification and State Space Representation 

If we consider only a uniaxial state of stress, the deviatoric applied and internal stress tensors become 
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Sii - 3 


and 


a,; = 3 


2<r 0 0 

0 — er 0 

0 0 -tr 


2a 0 0 

0 —a 0 
0 0 -a 


such that 


F = (^- v ) 


(35) 

(36) 


G = 


cr 


and the flow and evolutionary laws, respectively, become 


e J =^F' 

K 


and 


, (<r - a) 

| o-- a | 



(37) 

2R a B 0 ^ 2 «+!\ 

if 

(a — a)a > 0 

(38) 


(a- 

a)a < 0 

(39) 


Utilizing the concept of the state space representation [10] through [12], we can examine and illustrate 
geometrically the various features of the present unified viscoplastic model and discuss its characteriza- 
tion. The pertinent state space is constructed by combining and rewriting equations (7), (37), and (38), 
in standard vector form as follows, and plotting the corresponding Q fields or relaxation trajectories. 
That is, 


$=G (<r,a)+Z)(e) 


or, for example, during loading, 

UK 


-E t 


n{ 1 - §Ca 2 ) [<* - 0 II i 1 II Hv[Y]a - (*) 2,+1 ] 


)*{'.*} 


where 


;/ ii_ t 


H pn 
K 


n = 


c = 


3B 0 (1 + B\p (-f ) 2p 2 ) 
/c 2 (l + B,p (f ) 2p-2 (6p — 5)) 


(40) 

(41) 

(42) 

(43) 
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Alpha (MPa) 

Figure 2: State space associated with the material parameters in Table 1. 


and e is the total uniaxial strain rate. These relaxation trajectories are obtained from equations (37) 
and (40) by imposing a constant total strain rate (i.e., e= 0), thereby implying that S=Q (<?, <*)• 

Figure 2 illustrates the state space representation associated with equations (37), (38) and (39), and 
the material parameters given in Table 1. Immediately, a number of qualitative features of the first 
and third quadrants of the state space can be deduced from equations (37) and (40). For example, the 
locus of steady state (or fixed) points in creep, denoted by line ab in Fig. 2 is obtained by imposing the 
condition S = 0 on equations (40). Mathematically, such a condition is obtained when 


<r = (1 — 0Hv\Y])a + k 


1 + 


2R a B 0 /a\2»+i' 

Hisgnie 1 )- 0Hv[Y]a) W 


J.T 

n 


(44) 


and is denoted as the locus of steady state creep. The associated steady state creep rate is 


2R~Bo Co\ 2 «+ 1 

t — £-* \zl 

(sgn(^)-pHv[Y)a) 

Secondly, the locus of maxima in stress (i.e., <r = 0) during relaxation (e = 0) occurs when 


(45) 


k l = -1 = 0 (46) 

E 

This is satisfied when 

<r = k + a(l - fiHv[Y]) (47) 

and is denoted by line ac in Fig. 2. The points on lines ac do not correspond to fixed points in the state 
space, as do those associated with line ab, since the internal stress continues to evolve, that is, 
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Constants 

Units 

TIMETAL 21S 

Viscoplastic 



K 

MPa 

5.86 

V 

MPa/sec 

5.52xl0 -9 

n 

- 

3.3 

Bo 

MPa 

5.86xl0 -4 

B i 

- 

0.05 

P 

- 

1.8 

q 

- 

1.35 

R a 

1/sec 

O.lxlO -5 

P w 

- 

0.01 

Elastic 



E 

MPa 

80,671 

V 


0.3 


Table 1: Isothermal material parameters for TIMETAL 21S at 650°C 



Note that the inelastic strain is zero along and within the zone denoted by line ac, as in this region 
F < 0. Finally, Fig. 2 illustrates clearly that the second and fourth quadrants as well as the lower 
portions of quadrants one and three contain distinctly different relaxation trajectories resulting from the 
imposition of the internal unloading criterion. 

4*2 Characterization 

The most important, and often times most difficult, aspect of modeling the behavior of a given material 
at elevated temperature is obtaining the required material functions, e.g., f(F) and g(G), and material 
parameters. The difficulty associated with this process typically stems from not only the variety in 
mathematical forms for the material functions (e.g., power law, exponential, hyperbolic sine, etc.), but 
also the fact that given the material functions there is not a unique set of material parameters for any 
given load path. Therefore, numerous iterations and difficult compromises are required before a final set 
of material parameters, for the assumed material functions, can be obtained. 

With the present isothermal GVIPS model, nine material parameters need to be determined, as given 
in Table 1: three are associated with the flow law (i.e., «:,/*, and n); three with the nonlinear hardening 
operator, (i.e., S 0 , Bi, and p); two with the internal thermal recovery term (i.e., R a and q)\ and one 
with the dynamic recovery term (i.e., /?). Both elastic and inelastic material parameters given in Table 1 
correspond to the reference temperature (T re y) of 650°C. This temperature was chosen as the reference 
temperature because at 650 a C the material exhibits cyclically neutral behavior and thus supports the 
use of a purely kinematic model. Here a brief outline of the characterization process is given, with a 
more complete discussion regarding both exploratory and characterization testing given by Arnold [10] 
and Castelli et al. [19]. All tests addressed are uniaxial isothermal experiments, thus implying that the 
multiaxial material constants are typically generalized from their uniaxial counterparts. This need for 
generalization is precisely why a consistent multiaxial theory, such as that developed from a potential 
formulation, is imperative. The available tests for use in characterizing the current model are, three 
relaxation tests, three creep tests, one multiple-step creep test, and three tensile curves at different total 
strain rates. 

Initially, the yield stress is unity and therefore the radius of the threshold surface is determined by 
the drag strength and represents that value of stress below which only elastic behavior is "observed , 
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i.e., no measurable inelastic behavior occurs. The actual numerical value employed will obviously be 
non-unique as it will be dependent upon the sensitivity of the experimental equipment and/or definition 
of inelasticity employed. A preferred approach to arrive at a value for is to conduct a sequence of 
creep or relaxation probing tests to determine the maximum value of stress for which no time dependent 
behavior (c 1 = 0) is observed. An alternative approach and the one taken here (given the availability of 
three relaxation test results) is to use the minimum stress value attained at the end of a 24 hour period 
as our maximum allowable drag strength (i.e., k < 7 MPa) since all relaxation trajectories in the state 
space tend toward the origin (i.e., a = k, a = 0, see Fig. 2). 

Given k and taking /? « 0 (because at T re / the thermal recovery mechanism dominates) the material 
parameters fi and n may be obtained from either variable strain rate tests or stress transient dip tests, 
as discussed in Arnold [10]. Here, given a multiple step creep test which is a type of variable strain rate 
test, data pairs can be obtained at the beginning and end of each step in stress, e.g., (^i,f{), (<^ 2 >^ 2 )> 
respectively. Assuming that the internal stress a is equal to a* both before and after the step in stress 
(i.e., points 1 and 2) the following relation may be obtained from eliminating a* from the corresponding 
flow laws at points 1 and 2 (see equation (36)); that is, 




(49) 


where the only unknowns are and n. Employing a sufficient number of data pairs and a nonlinear 
regression scheme, optimal values for and n can be found. 

Next, assuming that hardening dominates the early portion of a creep test (primary creep stage) the 
thermal recovery term in the evolutionary law may be neglected. Consequently, the internal stress rate 
of equation (38) can be expressed as 


a =lH(l -|caV ( 50 ) 

Z o 

given that j3 & 0. Furthermore, an expression for a in terms of e 1 , e 1 and a can be obtained from 
differentiation of the flow law in equation (37), that is, 


a — 



(51) 


Thus, equating equation (50) and (51), assuming that k,/j, and n are known, considering data pairs 
(<r, c 1 , e 7 , a) taken from primary creep data, and using again a nonlinear regression scheme, the 

material parameters Bo,B\ y and P can be determined. Given multiple creep curves a number of sets of 
B 0 , B\> and p may be obtained and then a weighted averaged can be computed to give a single optimum 
set. 

Finally, the material parameters associated with thermal recovery ( R a and q) may be determined 
using the data pairs (<r 5 , c{) associated with the constant strain rate portion of the creep curve, and the 
corresponding internal stress (that is calculated from the flow law) at this assumed steady state (o,), 



14* 


(s< 


1 - 0Hv[Y] 


(52) 


Upon obtaining these data pairs, one may take the natural logarithm of both sides of equation (45) and 
determine R a and q from the ordinate intercept and slope, respectively, of the resulting ln-ln plot. 

The above outlined procedure can easily be applied to obtained an initial set of required material 
parameters. This set must, however, be further optimized to obtain a final set of material parameters 
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Figure 3: GVIPS correlation with experimental tensile data [18] considering the effect of strain rate. 

that will allow accurate correlation of the model over a ’’wide range of loading”. In the present case, the 
authors are primarily interested in the short term response (< 7200s) and small strain 1%) regime of 
the titanium alloy TIMETAL 2 IS. The resulting optimized material parameter set for TIMETAL 2 IS is 
given in Table 1, 

The corresponding correlation for, 1) the tensile response of TIMETAL 21S over three orders of 
magnitude in total strain rate (i.e., e = 8.33xl0 -4 , 8.33xl0” 5 , 8.33xl0 -6 ) , 2) short term creep at three 
stress levels (i.e., a = 72, 110, and 128 MPa), and 3) relaxation at 238 MPa are shown respectively in 
Figs. 3, 4 and 5, wherein the symbols denote experimental data and the lines the model correlations. 
In the authors’ opinions the model does a very good job, given a single set of constants, for the wide 
variety of loading conditions examined. It is important to remind the reader that if one were to focus 
on only one type of loading condition, superior correlations could be achieved for that class of loading; 
however, predictions of other classes of loading would severely suffer. For example during the above 
characterization process excellent tensile correlation was obtained but at the expense of poor creep and 
relaxation behavior for a given set of material parameters. Conversely, when the model was calibrated 
for creep response poor tensile and relaxation behavior were similarly predicted. 

5 GVIPS Predictions 

The present GVIPS model, characterized to represent the behavior of TIMETAL 21S at 650°C, is 
now exercised and its predictive capability assessed. The assessment begins by considering relaxation 
tests performed at stress levels above and below that used in the characterization of the model. The 
corresponding short term (i.e., < 800s) and mid term (i.e., < 7200s) relaxation responses are shown in 
Figs. 5 and 6, respectively. Again, it is apparent that overall agreement with short term behavior is 
good while that of mid term is only reasonable. Note how the initial stress rates are accurate at the 
higher stress levels and yet diminish more rapidly than do the actual experimental values. 
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Figure 4: GVIPS correlation with experimental data of TIMETAL 21S under creep. 



Figure 5: Short term relaxation response of TIMETAL 21S: GVIPS simulation versus experimental 
results. 
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Figure 6: Mid term relaxation response of TIMETAL 21S: GVIPS simulation versus experimental results. 

Figure 7 illustrates the cyclic behavior under strain control given a total strain rate of 8.33 x 10 -5 , 
wherein the symbols represent the experimental results and the solid line the GVIPS prediction. Again, 
the agreement is adequate with approximately a 20% over-prediction of the maximum stress at ±0.01 
strain. Although not shown in Fig. 7, note that two cycles are required with the GVIPS model before 
the hysteresis loop will close. 

A multiple step creep test is considered next, where the resulting inelastic strain versus time (at 
constant stress) is shown in figure 8. Here the GVIPS (solid line) simulation is seen to significantly 
under-predict, the accumulated inelastic strain as compared with the experimental data (denoted by the 
symbols). This is not surprising, however, given the Fact that the model under-predicts the creep behavior 
at 72 MPa (see Fig. 4) and that this test was performed under constant load and not constant stress 
as was the simulation. As a consequence of the constant load controlled creep test, the experimental 
results have both additional geometric effects and possibly creep damage (thereby leading to tertiary 
creep) included in them. Neither of these factors have been included in the numerical simulation. The 
influence of these factors is particularly evident during the third step, where the inelastic strain rate 
during the experiment is significantly greater than that simulated (cf. <r = 128 MPa in Fig. 4). 

Lastly, a classic plasticity-creep interaction experiment was performed and simulated to illustrate 
the need for a unified viscoplastic formulation. Here the material is subjected to an overload prior to 
performing a creep test at a lower stress value, as shown in Fig. 9. Typically, as is found here, the 
initial creep response following the overload is significantly reduced when compared to that produced 
from a pure creep test at this lower stress amplitude, see Fig. 10. Once again comparing the numerical 
simulation (of the GVIPS model) to that of the experimental results one observes reasonable overall 
agreement. Although, the amount of total accumulated creep strain at the end of the test, i.e., 72000 
seconds, is significantly under estimated, the initial reduction in primary creep rate, due to the overload, 
is extremely accurate during the first 800 seconds and good over the first hour of creep loading, thus 
fulfilling our primary objective. 
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Figure 8: Multistep creep test: GVIPS and BP simulations versus experimental data 
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Figure 9: Plasticity - Creep interation test: stress-strain reponse simulation versus experimental data. 
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Figure 10: Plasticity - creep interation test: illustrating the influence of prior history on inelastic re- 
sponse. 

6 GVIPS vs. BP Model 

Reviewing figures 3 through 10 it is evident that the present GVIPS model in conjunction with the 
material parameters of Table 1 does a very good overall job of simulating the short term and small strain 
behavior of TIMETAL 21S at 650°C. In an attempt to assess the relative quality of these correlations 
and predictions as well as the numerical efficiency of the present model, an alternative and commonly 
employed unified viscoplastic model was implemented into our in-house Inelastic Deformation Analysis 
Code (IDAC). The full multiaxial, nonisothermal formulation of the Bodner-Partom (BP) model, given 
in [5] was chosen as the alternative model because it too has been characterized for the titanium matrix 
TIMETAL 21S. The specific tests employed to characterize this model are unknown to the authors, thus 
all of the previous loading conditions will be investigated and assumed to be predictions. 

Figures 11 through 15 illustrate the tensile, creep, cyclic, an <T relaxation behavior predictions of the 
BP model at 650°C. The corresponding GVIPS simulations are also included in these figures for ease of 
comparison. Examining Fig. 11 it is apparent that at the higher strain rates (8.33xl0~ 5 and 8.33xl0~ 4 ) 
the BP model does a better job of predicting the saturated stress values at 0.01 strain than does the 
GVIPS model. However, at the lower strain rates the BP model does a poorer job in representing the 
knee (or transient portion) of the tensile curves. This poorer simulation by the BP model of the time 
dependent behavior is further illustrated in Figs. 12, 14, and 15, wherein the creep and relaxation 
behavior is compared with the GVIPS simulations. Clearly, both the primary and secondary creep rates 
differ greatly and disagree with the experimental observations. Similarly, the short and mid term BP 
relaxation response is inaccurately predicted in Figs. 14 and 15 as compared with the more accurate 
GVIPS simulations. This inferior performance in creep and relaxation is not unexpected as the BP 
model is known for its ability to predict the rate-dependent tensile and cyclic response of a material 
and not the time-dependent creep and relaxation behavior. The significance of this failure to adequately 
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Figure 11: GVIPS and BP simulation of the influence of strain rate on the tensile response. 


capture the time-dependent behavior is further emphasized in both the multi-step creep test shown in 
Fig. 8 and the plasticity-creep interaction test shown in Figs. 9 and 16 wherein the BP model poorly 
represents the primary creep (transient or short term response) following a change in the applied loading. 

The ability of the BP model to simulate the rate-dependent cyclic behavior is shown in Figs. 7 and 
13, where both the BP and GVIPS responses are in reasonably good agreement, except for the predicted 
maximum stress amplitudes achieved, where they differ between 10% and 30% of one another. Also, 
note that in Fig. 13 a numerical problem was observed in the BP model during unloading at the fastest 
strain rate of 8.33x10 “ 4 . No special numerical techniques where implemented to attempt to overcome 
this difficulty as suggested in the work of Kroupa [5], because a strict comparison in the numerical 
integration of these two models was desired. 

Finally, a comparison of the numerical efficiency of both models is made, relative to the total number 
of integration steps required to complete the above simulations. This comparison is completely legitimate 
as both multiaxial models were implemented into the same computer code, integrated using the same 
self-adapting fourth order Runge-Kutta and Adams Bashforth predictor-corrector integration algorithm, 
and executed on the same Sun IPX workstation. Table 2 contains the number of total integration steps 
required to perform each loading history discussed. Clearly, the GVIPS model is extremely easy to 
integrate (due to the use of power law material functions as opposed to exponential ones and the 
potential structure itself) when compared with the BP model and therefore would provide significant 
computational savings in a large-scale structural finite element analysis. Note that in the case of tensile 
and cyclic behavior both models give similar accuracy and yet the GVIPS formulation provides at least a 
factor of 2 improvement in numerical efficiency over the BP model. Furthermore, in the case of creep and 
relaxation where the simulation accuracy of GVIPS is clearly superior to the BP model, the numerical 
efficiency is also superior, particularly in the case of creep where numerical stiffness can become a problem 
[21]. Thus the GVIPS (generalized viscoplasticity with potential structure) formulation provides not 
only accurate predictive capability, but also superior numerical performance over other non-potential 
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Figure 15: GVIPS and BP simulation of mid term relaxation response. 


21 






Creep Time (s) 


Figure 16: GVIPS and BP simulation of inelastic (creep) response following overload shown in Fig. 9. 


Type of Test 

Integration Steps 

Speed-up 


GVIPS 

BP 


Tensile (1/sec) 
e = 8.33xl0 -4 

170 

338 

1.99 

= 8.33xl0 -5 

196 

619 

3.16 

= 8.33x10"® 

231 

992 

4.29 

Cyclic (1/sec) 
e = 8.33xl0 -4 

761 

15007 

19.72 

= 8.33x10"® 

828 

2727 

3.29 

= 8.33x10-® 

938 

4937 

5.26 

Creep (MPa) 
a = 72.4 


758 

2.26 

= 110 



34 

= 128 

364 

32621 

89.6 

Relaxation (MPa) 
<j = 345 

372 

624 

1.68 

= 238 

365 

519 

1.42 

= 103 

358 

382 

1.07 

4 Step Creep Test 

1189 

6631 

5.57 

Plasticity /Creep 

476 

8839 

18.5 


Table 2: Numerical Comparison of the GVIPS and BP Models 


22 
















formulations, the BP model being an example. 


7 Summary 

A fully associative, multiaxial, isothermal, nonlinear kinematic hardening viscoplastic model has been 
presented. It contains in general, three internal state variables (two scalars and one tensor) and both 
thermal and strain- induced recovery mechanisms. The two, non-evolving, scalar internal state variables 
are associated with dislocation density and are defined as the drag and yield stress. The evolving 
tensorial variable known as the internal (or back) stress is a second order, traceless, symmetric tensor 
and is associated with dislocation substructure. A unique aspect of the present model is the inclusion of 
non-linear hardening through the use of a compliance operator, derived from the Gibb’s potential, in the 
evolution law for the back stress. This non-linear tensorial operator is significant in that it allows both the 
flow and evolutionary laws to be fully associative (and therefore easily integrated) and greatly influences 
the multiaxial response under non-proportional loading paths. In addition to this nonlinear compliance 
operator, a new consistent, potential preserving, internal strain unloading criterion has been introduced 
to prevent abnormalities in the predicted stress-strain curves, which are present with nonlinear hardening 
formulations, during unloading and reversed loading of the external variables. The resulting isothermal 
unified viscoplastic model was then characterized for the titanium based alloy TIMETAL 21S at 650°C. 
Results illustrated the good overall correlation and predictive capabilities of the model for a wide range 
of loading conditions, that is tensile, cyclic, creep, and relaxation tests. The proposed model was then 
compared with a commonly accepted and employed version of the Bodner-Partom viscoplastic model 
and was found to be superior both in its predictive capabilities and numerical performance. 
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